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ABSTRACT 

Although the influence of magnetic fields is regarded as vital in the star formation 
process, only a few magnetohydrodynamics (MHD) simulations have been performed 
on this subject within the smoothed particle hydrodynamics (SPH) method. This is 
largely due to the unsatisfactory treatment of non-vanishing divergence of the mag- 
netic field. Recently smoothed particle magnetohydrodynamics (SPMHD) simulations 
based on Euler potentials have proven to be successful in treating MHD collapse and 
fragmentation problems, however these methods are known to have some intrinsical 
difficulties. We have performed SPMHD simulations based on a traditional approach 
evolving the magnetic field itself using the induction equation. To account for the nu- 
merical divergence, we have chosen an approach that subtracts the effects of numerical 
divergence from the force equation, and additionally we employ artificial magnetic dis- 
sipation as a regularization scheme. We apply this realization of SPMHD to a widely 
known setup, a variation of the 'Boss & Bodenheimer standard isothermal test case', 
to study the impact of the magnetic fields on collapse and fragmentation. In our sim- 
ulations, we concentrate on setups, where the initial magnetic field is parallel to the 
rotation axis. We examine different field strengths and compare our results to other 
findings reported in the literature. We are able to confirm specific results found else- 
where, namely the delayed onset of star formation for strong fields, accompanied by 
the tendency to form only single stars. We also find that the 'magnetic cushioning ef- 
fect', where the magnetic field is wound up to form a 'cushion' between the binary, aids 
binary fragmentation in a case, where previously only formation of a single protostar 
was expected. 

Key words: magnetic fields - MHD - stars: formation - ISM: clouds - ISM: magnetic 
fields 



1 INTRODUCTION 

Magnetic fields, besides self- gravity, radiation and turbu- 
lence, are usually regarded as being the most fundamental 
constituents needed to describe star formation (for recent 
reviews of theoretical aspects see, e. g., |Mac Low & Klessen 
2004||McKee & Ostriker|2007[ and references therein). Espe- 
cially the effects caused by magnetic fields, supported by in- 
creasing observational evidence of magnetic field structures 
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in the interstellar medium (ISM) and molecular clouds (see, 
e. g., Heiles &; Crutcher|2005[ and references therein), came 
into the focus of interest within the past decade. 

Numerical simulations, however, were not able to han- 
dle the full complexity of the physical processes connected 
with star formation for a long time, and were therefore lim- 
ited to pure hydrodynamical, self-gravitating investigations. 
But in recent years the situation changed dramatically. Eu- 
lerian codes, which were always being able to handle mag- 
netohydrodynamics (MHD) with good accuracy, became, 
with the advent of adaptive mesh refinement (AMR, see 
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Berg er &; Colella|r 989 ) , able to handle protostellar collapse 
with MHD. Examples include the investigations by Ziegler 
(2005} and |Fromang et aL (2006), who included collapse 
problems based on variations of the standard isothermal 
test case ( |Boss &; Bodenheimer 1979), to test their AMR 
codes, NIRVANA and RAMSES ( |Teyssier|2002[ ), respectively. 
They found a strong influence of magnetic fields on proto- 
stellar fragmentation in the limit of ideal MHD, that is, in 
a medium with infinite conductivity. Further investigations 
were performed by Commergon et al. (2010), emphasizing 



the importance of considering the combined effects of MHD 
and radiative transfer on collapse and fragmentation. 



Al so based on AMR, |Machida et al.| ( |2QQ4| ) and |Machida 
et al. ( 2005a|b ) performed several collapse simulations in 



ideal MHD, and found that fragmentation was suppressed 
by the magnetic field, but occurred still. More recently, these 
studies were extended to the formation of metal-free Popu- 



In a series of publications, 


Hennebelle & Fromang 


(2008) and Hennebelle & Teyssier 


(2008) also investigated 



the effect of magnetic fields on the collapse of dense molecu- 
lar cloud cores. The former concentrates on magnetic brak- 
ing and launching of outflows, where models with differ- 
ent magnetic field strengths are considered. For weak fields, 
they found negligible magnetic braking and thus formation 
of a centrifugally supported disc which in turn triggers a 
slowly expanding magnetic tower. For higher magnetic field 
strengths, they did not find formation of a centrifugally sup- 
ported disc as a consequence of strong magnetic braking and 
collapse along the field lines. The latter publication however, 
focuses on fragmentation where a perturbation is added to 
the same setup as in the former work. Here, at weak field 
strengths, the centrifugally supported disc, which fragments 
in the hydrodynamic case, is found to remain stable and ax- 
isymmetric. For strong magnetic fields, again, no centrifu- 
gally supported disc is found because of magnetic braking 
and fragmentation is only found for strong initial perturba- 
tion amplitudes. 

Jets and outflows, however, are closely associated 
with star formation, and so many MHD simulations have 
been performed to investigate these phenomena, which are 
thought to be driven by coupling to magnetic fields. The 
study of jets and outflows is a field on its own right, so we 
and references therein, for a corn- 



refer to |Banerjee| ([2009} 
prehensive discussion. 

For more than 20 years, there has been an attempt to in- 
clude MHD in smoothed particles hydrodynamics (SPH) for 
use on collapse problems, starting with the work of Phillips 
(1986a b). However, his code lacks important algorithmic 
features developed afterwards and regarded as vital ingredi- 
ents in SPH codes today, like adaptive smoothing lengths. 
More seriously, he considered only non-rotating clouds, so 
the results, which show no fragmentation either in the mag- 
netized nor in the non-magnetized clouds, have to be taken 
with a grain of salt. 



Hosking & Whitworth ( 2004 ) used a different approach 
based on a two- fluid formalism. This allowed non-ideal ef- 
fects, such as ambipolar diffusion, to be taken into account, 
and thus enabled them to start from a subcritical rotating 
cloud core which, after following the evolution for some time, 
turned supercritical as result of diffusion. Their general con- 
clusion was, that magnetic fields inhibit fragmentation. But 



their implementation suffered from non-zero divergence of 
the magnetic field, which did not allow them to follow the 
evolution for a long time. 

An important step forward was the study by |Price fe| 



Bate (2007). Their implementation is based on Euler po- 
tentials ( |Stern|1970| [Rosswog fe Price|2007| ), which are free 



of physical divergence by construction. In their work, they 
considered two well known models, namely the axisymmet- 
ric collapse of a homogeneous density sphere and a vari- 
ant of the Boss & Bodenheimer (1979| 'standard isothermal 
test case' with initial m — 2 perturbations in density. They 
found, that stronger magnetic fields caused delays to the 
collapse, since the additional magnetic pressure provides ad- 
ditional support against gravity. Furthermore, they pointed 
out, that potentially crucial effects on discs might be caused 
from this delay, since the rate of mass infall onto the disc is 
reduced in this case. With respect to the perturbed clouds, 
these authors drew the main conclusions that magnetic fields 
might not be a serious problem to binary formation but that 
they suppress fragmentation. The latter is, contrary to pre- 
vious results reported in the literature, attributed to the 
additional support by magnetic pressure, rather than mag- 
netic tension forces or magnetic braking. 

In their following works, they turned to magnetic fields 
in cluster formation ( Price &; Bate|2008 ). Using a barotropic 
equations of state, they found differences compared to pure 
hydrodynamical runs. Especially a significant influence to 
the star formation rate was reported by these authors. They 
performed further investigations by replacing the equation 
of state with a radiation transfer treatment, based on the 
flux-limited diffusion approximation (Price feBate|2009) . As 
a main result of their investigations, the authors conclude 
that the net result of magnetic fields and radiative transfer 
is able to explain the inefficiency of star formation, with a 
star formation rate of 10 per cent per free-fall time, which 
is in good agreement with observations. 

However, despite their obvious success, it must be noted 
that Euler potentials have limitations on their own. The 
most serious one is that magnetic helicity is constrained to 
be zero, since the vector potential is always perpendicular 
to the magnetic field. In practice, this means that certain 
field configurations, namely such that are mult i- valued as, 
e. g., combinations of poloidal and toroidal geometries, can- 
not be represented by Euler potentials. From this it follows, 
that such configurations also cannot be generated during 
a simulation, making it impossible to study certain physi- 
cal processes, e. g., winding up of magnetic fields as found 
in dynamos or in protostellar outflow phenomena. Another 
limitation was pointed out recently by |Brandenburg| ( |2010| ), 
who stressed the fact that Euler potentials are not able to 
deal with non-ideal MHD since even a small amount of dif- 
fusivity prevents convergence to the correct solution. 

So other ways of dealing with the magnetic fields need 
to be considered, and the research on this topic is still on- 
going. Another promising idea is an approach based on the 
vector potential, and indeed for spacial dimensions smaller 
than three, good results have been obtained by Price (2 010| ). 
This is due to the fact that in these cases, the vector poten- 
tial is in fact mathematically equivalent to a formulation 
with Euler potentials. However, |Price (2010) also showed 



that the vector potential formulation was not even able to 
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the development version of GADGET-3. We used the code 
only in non-expanding, Newtonian space, so all equations 
referring to implementation details are lacking cosmological 
parameters and extensions. 

In problems related to star formation, physical quanti- 
ties vary over several orders of magnitude. Therefore, spatial 
and temporal adaptivity must be guaranteed within the sim- 
ulation. The former is done using individual and adaptive 
smoothing lengths, where for each particle i the equation 



4?r 3 



Nrm 



(1) 



is solved iteratively with the density. Here, hi is the parti- 
cles smoothing length, N the number of neighbours, m, is 
the particle mass and pi is the density which is calculated 
according to 



Figure 1. The running average for the mean density at time 
t = is shown, where every cross indicates the mean within a 
bin of length R/n^i n and the corresponding standard deviation 
is shown as error bar. We chose ribin = 50 for this analysis. The 
red line is a fit of the points to a constant function, obtaining a 
value of 0.998 and thus a deviation from the analytical po by only 
0.2%. 



handle standard test cases in three dimensions, causing him 
to suggest not to use this approach in an SPH context. 

In this paper, we follow a different approach using the 
MHD implementation into the widely used GADGET code 
(Springel et al. 2001| ||Springel 2005), which has been ap- 
plied successfully to several problems in galactic astrophysics 
( Dolag fe Stasyszyn|[2009l |Kotarba et"aT||2009[ |2010) . We 
use the traditional approach for ideal MHD using the in- 
duction equation, but subtract the non- vanishing divergence 
term from the force equation to ensure numerical stability 
(B0rve et al. 2001). As regularization scheme, we employ 



time-dependent artificial resistivity, introduced by |Price~fe] 



Monaghan (2005). We show, that this method produces ac- 
curate results which are not corrupted by non- vanishing di- 
vergence and compare well to some of the findings of |Price| 



& Bate (2007). However, at higher field strengths, we see 



noticeable deviations, the most prominent being the 'mag- 
netic cushioning effect', where the magnetic field is wound 
up due to the rotation of the cloud and forms a cushion be- 
tween the protostars. The latter effect thus aids binary star 
formation in the case of a mass-to-flux ratio of 4 (in critical 
units) where Price & Bate (2007) had found just a single 
protostar. 



2 METHOD 
2.1 Code 

For the (magneto-) hydrodynamical simulations presented in 
this work, we use the GADGET code ( jSpringel et al.||2001| 
Springel 2005), a tree-based, massive parallel code utilizing 
the SPH method (for recent reviews of SPH see, e. g., |Ross-| 
|wog||2009 Springel 201 0|). The simula tion results presented 
here were, as in |Dolag &; Stasyszyn (2009), obtained with 



Pi y^mjW(rjj,hi) 



(2) 



where W is the cubic spline kernel (Monaghan & Lattanzio 
|1985| . The dynamical equation 



dvi 



(hyd) 



3 = 1 



Pi 



>5> 



(3) 



has been derived using a variational principle (e. g. | Springel] 
|&; Hernq"u ist 2002), and so the so called "grad h" correction 
terms 



fi 



1 + 



hj dpi 
Spi dhi 



(4) 



which account for the derivative of the kernel with respect to 
the smoothing length, are included by construction, ensuring 
energy and entropy conservation to time step accuracy. 

To allow for accurate shock capturing, artificial viscos- 
ity is needed. The contribution of the viscous term to the 
particle acceleration is given by 



dvi 
~dt 



(vise) 



-^rrijUijViWij, (5) 



where Wij = (Wi + Wj)/2, and is the viscous tensor 
which is defined as 



2/V 



o, 



for Vi 



for Vi 



<0, 



(6) 



> 0. 



Vj and Vij 



where p^ = (pi + pj)/2, Vij = 
This term was derived in close analogy to Riemann solvers 
(see Monaghan 1997) and includes the signal velocity 



V ii =Ci+ Cj 



f3vij 

dPi/dpi = 



(7) 



with the sound speed d 

(3 we use, as suggested by Dolag & Stasyszyn (2009), the 



^jP^JjH . For a and 



values 2.0 and 1.5, respectively. Additionally, we would like 
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Figure 2. Simulation results without magnetic field, i. e. the pure hydro dynamical case. In horizontal direction, x-coordinates are 
plotted, and y-coordinates are plotted vertically. Shown is the column density, here integrated along the z-axis, in physical coordinates. 
The panels display a sequence of increasing time, measured in in units of free-fall times, % = 2.4 x 10 4 yr. 



to mention that we did not use the viscosity limiter intro- 
duced by Balsara (1998), since it is very likely responsible 



for introducing numerical artefacts in magnetic field growth, 
as observed in other work ( Kotarba et al.||2009 |2010| ). 

It should be noted, that artificial viscosity is a source 
of entropy (so is artificial resistivity, see below), which is 
generated at a rate dAi/dt, where A = P/p 1 is the entropic 
function. Since we use a barotropic equation of state, thus 
calculating the pressure directly as a function of density, an 
explicit consideration of the entropy production is not nec- 
essary and so we do not make use of the entropy treatment 
in GADGET. For a detailed discussion of the entropy for- 
mulation in SPH and implementation details, we refer to 
Springel fe Hernquist| ( |2002| ) and |Springe"I| ( |2QQ5| ). 

The timestepping scheme is adaptive and an individual 
timestep for each SPH particle is chosen as a minimum of 
two criteria, 



Ati — min 



2r\Ei Co 



thi 



(8) 



The first of these is based on a particles acceleration where r\ 
is an accuracy parameter, with a numerical value of 10 -3 in 
this work, and e the gravitational softening length, while the 
second is a Courant-like criterion that is needed to ensure 



numerical stability. For the sink particles (see section 2.5), 
only the first criterion is used. 



2.2 Magnetohydrodynamics 

For the work presented here, we use the SPMHD implemen- 



tation into GADGET described in detail in Dola g &; Sta-| 
|syszyn| (2009). In the latter work, extensive tests on the 
reliability of the algorithms have been performed, among 
them some well known standard test cases typically used in 
the literature. These include the shocktubes considered in 
the work of Ryu &; Jones| ( |1995| ) , the fast rotor by|Balsara 
& Spicer (1999), the Orszag-Tang vortex (lOrszag & Tang 



1979 ) and a variation of the strong blast test (e. g. |Balsara 
fe Spicer|1999l ). It was shown, that this implementation per- 



forms very well in general. Here we repeat the fundamental 
parts of the implementation and refer to |Dolag &; S tasyszyn 
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Figure 3. Simulation results with magnetic field parallel to the rotation (z) axis, with an initial field strength of Bo = 40.7 \xG. This 
corresponds to a mass-to-flux ratio of M/<l> = 20, in multiples of the critical value. As previously, in horizontal direction, ^-coordinates are 
plotted, and y-coordinates are plotted vertically. Shown is the column density, here integrated along the z-axis, in physical coordinates. 
The panels display a sequence of increasing time, measured in in units of free-fall times, % = 2.4 x 10 4 yr. 



( 2009 ) for a more detailed discussion on algorithms and per- 
formance in test cases. 

For the correct capturing of shocks in the hydromag- 
netic case, it is essential to assign the correct artificial vis- 
cosity to the particles. Therefore, the sound speed a in eq. 
([7| is replaced by the speed of the fastest magneto-sonic 



(B) 

M V J 



V2 



(c? + t&) + */(c?+*£) 2 -4 



Mo Pi 



1/2 



(9) 



which enters also in the timestep criterion via the Courant 
condition. 

The contribution of the magnetic field to the accelera- 
tion is given by 



dv k \ (B) ldM kl 



dt 



p dx l 



(10) 



where M is the magnetic stress tensor (Phillips &; Mon 



agha n|1985 ) defined as 



Mo 



B k B l 



-\B\ 2 5 kl 
2 1 1 



(11) 



where, as in equation (10), the upper indices denote coor- 



dinates. A straightforward discretization of the equation of 
motion, which also can be derived from a Langrangian using 
a variational principle (Price & Monaghan 2004b), is 



dvj 
dt 



(B) 



Mo 



(12) 

and Vi, and the corresponding formulation for 
are abbreviations for M kl and Vf, respectively. 
While this formulation of the magnetic force conserves mo- 
mentum exactly, it is also known to be unstable to negative 



where M.% 
particle j, 



stresses causing the particles to clump ( Phillips & Monaghan 



6 F. Biirzle et al. 



1985). While many possible methods have been proposed in 



the literature to correct for this instability, most of them 
are rather impracticable or only of limited use, see |Dolag| 
& Stasyszyn (2009) for a detailed discussion. In this work, 



we choose the formulation introduced by B0rve et al. (2001 ) 



which subtracts the effect of any numerically non- vanishing 
divergence of the magnetic field. This is done by subtracting 



dvi 



Mo 



7=1 



P? 



•v,ir. • /;- ' -v.ir, 

(13) 



Dolag & Stasyszyn ( |2009| ) where this method was used 



throughout and gave excellent results. Furthermore, the ef- 
fects due to violation of momentum conservation have been 
shown to be negligible. 

In ideal MHD, that is, in a medium with infinite elec- 
tric conductivity, the magnetic field is advanced using the 
induction equation (Price & Monaghan 2004b) 



(IB 

dt 



= (B- V)v-B(V-v), 



(14) 



and its SPH discretization is given by 



dBj 
dt 



^rrij [Bi{vij • ViWi) - Vij(Bi • ViWi 



(15) 

where the "grad /i" correction terms are included for con- 
sistency (but note that this can not be derived from first 
principles) . 

Since a further important source of errors is the noise 
introduced by numerical fluctuations of the magnetic field, 
which originate in integration errors, a regularization proce- 
dure is required in the numerical scheme. We use artificial 
magnetic dissipation to regularize the underlying magnetic 
field. This is done in close analogy to the artificial viscosity 
by introducing a parameter a b that controls the strength of 



the dissipative effect. As in Dolag & Stasyszyn (2009) and 
Price & Monaghan ( 2005| , the dissipative term is included 
into the induction equation 



dBj 
dt 



(diss) 



2 ^ p.. 

7 = 1 H w 



(Bi — Bj)rij • V^W* 



(16) 

Note that as can be a constant or a time dependent quan- 
tity. In the latter case, a decay equation 



da B 
~df 



a B - a£ m 



(17) 



is evolved for each particle, where the source term S is given 
by 



So 



max(|V-B|,|V x B\) . 



(18) 



hi 



C max 



(19) 



where C is typically chosen in the same range as for the 
Courant timestep condition. In this work, we used the time- 
dependent version throughout, enforcing a maximum value 
of a# ax = 1.0 for the resistivity parameter. The value of 
the magnetic field constant /io was chosen such, that the 
magnitudes of the magnetic fields are given in Gauss. 



2.3 Comments on artificial viscosity and artificial 
resistivity 

While artificial viscosity is needed to allow for correct shock 
capturing and to avoid unphysical effects such as particle in- 
terpenetrations, this approach certainly has its well known 
weaknesses. One of the most prominent, pointed out (among 
others) by Agertz et al. ( 2007 ) , is the fact that with viscos- 



ity present in the system, one effectively needs to solve the 
Navier- Stokes and not the Euler equations. To avoid errors 
resulting from this, one method is the utilization of time 



dependent artificial viscosity ([Morris & Monaghan 1997) 



combined with a switch, so that viscosity can be reduced 
to a minimum if no sources of viscosity are present. How- 
ever, this approach is problematic when applied to collapse 
problems, since the usual switch is based on the condition 
V • v < 0, indicating the presence of a shock. In context 
of a self-gravitating collapse, this condition is also a sign 
for a convergent flow and thus the switch might erroneously 
respond in this case. So since we assume here, that the er- 
rors due to the latter effect would be more serious in our 
application than errors from not solving the correct Euler 
equations, we decided to use a constant artificial viscosity 
throughout. However, to reduce the effect of intrinsical nu- 
merical diffusion, we followed the conclusions drawn by the 



work of Attwood et al. (2007), and restricted the allowed 



range of the nearest neighbours N to a value smaller than 
one. In this work, we used N = 64 zb 0.3 throughout. 

The situation is similar with artificial resistivity. As al- 
ready mentioned, artificial magnetic dissipation is needed 
to deal with noise related to magnetic fields. However, this 
means that in principle we need to solve the equations of 
non-ideal MHD to account for the additional diffusivity 
added by artificial resistivity, so using the time-dependent 
formulation is likely to improve the situation considerably. 
Furthermore, a constant dissipation of, say, as — 1.0, which 
is needed in later stages of the collapse would introduce a 
large scale initial diffusion especially in the low density re- 
gion of our system and thus lead to a significant disturbance 
at early stages. 



2.4 Thermodynamics 

In our models, we use a piecewise equation of state, defined 
by 



A natural choice for the characteristic decay time-scale r 
is provided by considering the time a shock needs to travel 
through one kernel length and can therefore be written as 
( Price fe Monaghan||2004a ) 



Kp 1 . 



(20) 



Since this equation of state is barotropic, i. e. the pressure 
is a function of the density only, the energy equation needs 
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Figure 4. Simulation results with magnetic field parallel to the rotation (z) axis, with an initial field strength of Bo = 81.3 U.G. This 
corresponds to a mass-to-flux ratio of M/<& = 10, in multiples of the critical value. As previously, in horizontal direction, ^-coordinates are 
plotted, and y-coordinates are plotted vertically. Shown is the column density, here integrated along the z-axis, in physical coordinates. 
The panels display a sequence of increasing time, measured in in units of free-fall times, % = 2.4 x 10 4 yr. 



not to be solved explicitly. In this work, the adiabatic index 
is given by 



1, for p < yOcrit, 

7/5, forp>p C rit. 



(21) 



with p c 



10 14 gcm 3 . So for low densities the equation 



of state is isothermal and K = c s , while for high densities 
it is adiabatic assuming a diatomic gas with five degrees of 
freedom. In the latter case, K is chosen such that the pres- 
sure is continuous at the critical density, i. e. K — Cgp"^/ 5 . 



2.5 Sink particles 

In those high density regions that are going to form a pro- 
tost ar, particles also gain large accelerations which in turn 
leads to assignment of very small time steps to a small frac- 
tion of the particles present in the whole system. There- 
fore, the timestep is becoming prohibitively small, and ef- 
fectively causes the simulation to stall. Sink particles, first 



introduced by Bate et al. (1995), provide a way of solving 



this problem. When a certain threshold density is reached 
within some small region of space, characterized by the sink 



radius, the gas particles within this region are replaced by a 
non-gaseous particle that carries their masses and momenta. 
This particle interacts with other particles via gravity only, 
and is able to accrete further particles that cross its outer 
boundary. However, while making further evolution of the 
collapse accessible, the method comes with the burden that 
all information within the sink particle is lost. But after all, 
sink particles have proven to be a very useful subgrid model 
which has been successfully applied in many star formation 
related studies. 



Since the pioneering work of Krumholz et al. (2004), 



also Eulerian codes can benefit from sink particles, and re- 



cently Federrath et al. (2010) have accomplished a imple- 



mentation into the widely used FLASH grid code (Fryxell 
et al. 2000| . While the first implementation of sink parti- 
cles into GADGET-2 was done by |Jappsen et al. (2005), the 
FLASH implementation served as a prototype for our current 
implementation of sink particles in GADGET-3. 

In the studies presented here, we insert a sink particle, 
once a threshold density of p s = 10 -10 gcm -3 within an 
accretion radius of ~ 13 au has been reached. We chose the 
creation density to be far in the adiabatic regime, to ensure 
that sinks are only formed in regions where the collapse has 
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Figure 5. Simulation results with magnetic field parallel to the rotation (z) axis, with an initial field strength of Bo = 108.5 fiG. This 
corresponds to a mass-to-flux ratio of M/<J> = 7.5, in multiples of the critical value. As previously, in horizontal direction, ^-coordinates 
are plotted, and y-coordinates are plotted vertically. Shown is the column density, here integrated along the z-axis, in physical coordinates. 
The panels display a sequence of increasing time, measured in in units of free-fall times, % = 2.4 x 10 4 yr. 



advanced already several orders of magnitude thus ruling 
out artefacts by spurious sink formation. The gas that is 
approaching a sink particle later, is accreted if it is bound 
to it and further criteria are met, as described in detail in 



Federrath et al. (2010) 



The treatment of the magnetic field in sink particles, 
however, has the same limitations as pointed out in |Price ~&\ 
Bate (2007 2008), namely that magnetic field carried by ac- 
creted particles is discarded, so sink particles can not provide 
magnetic field driven feed back on the surrounding cloud. 



3 INITIAL CONDITIONS 



For comparison with Price & Bate (2007), we chose the same 
initial setup as in their work. The initial cloud core has a 
spherical shape with a radius R = 4 x 10 16 cm and a mass 
M = IM0, i. e. we adopt a constant initial density po = 
7 A3 x 10 _18 gcm -3 and a free-fall time 



tff = 



3tt 
32Gp 



: 2.4 x 10 yr. 



(22) 



The initial setup is realized by distributing the particles 
on a closed-packed lattice. This kind of particle distribution 
ensures very good settling properties with a very low ini- 
tial scatter. In Fig. |l]) we show the running average of the 
density, as obtained from an initial snapshot of the system. 
An analysis by fitting the points shows, that the deviation 
from the analytical density po over the whole cloud is only 
0.2%. So we conclude, that the influence of Poission noise, 



which, as discussed in Cartwright et al. (2009), imposes se- 
rious problems for SPH estimators, is effectively reduced in 
our initial conditions. Additionally, we would like to point 
out that the particle distribution is identical for every con- 
sidered model, as are all other physical parameters not re- 
lated to magnetic fields. 

The number of particles in the cloud are 300 914. As 
pointed out in Price & Bate (2007), these are about ten 
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Figure 6. Simulation results with magnetic field parallel to the rotation (z) axis, with an initial field strength of Bo = 203 |U.G. This 
corresponds to a mass-to-flux ratio of M/<1> = 4, in multiples of the critical value. As previously, in horizontal direction, ^-coordinates are 
plotted, and y-coordinates are plotted vertically. Shown is the column density, here integrated along the z-axis, in physical coordinates. 
The panels display a sequence of increasing time, measured in in units of free-fall times, % = 2.4 x 10 4 yr. 



times as many particles as required from the Jeans resolu- 
tion criterion (Bate fc Burkert|1997 ) for the chosen equation 
of state. For the MHD calculations, this cloud was embed- 
ded into a uniform, low-density medium with a temperature 
30 times higher than within the cloud, so that cloud and 
medium are initially in pressure equilibrium. This approach 
has the advantage, that the magnetic field lines behave reg- 
ularly at the cloud boundaries and thus cloud particles are 
prevented from being ejected into space by magnetic forces, 
which otherwise would be induced by ill-defined behaviour of 
the magnetic field at the cloud surface. The ambient medium 
is represented by 146 074 particles which does not add sig- 
nificant extra computational cost. Note, however, that for 
larger systems typically studied in the context of star clus- 
ter formation, this approach is not efficient any more, and 
another strategy has to be employed by, e. g., removing par- 
ticles far away from the region of interest. To avoid dilution 
of the medium, the whole system is placed into a cubic box 
with periodic boundary conditions and a side length of four 
times the cloud radius. 



We consider a variation of the 'standard isothermal test 

( 19791), but in an SPH con- 



Boss & Bodenheimer 



Bate Burkert (1997)] which, in the hydrody- 



case [see 
text also 

namical case, is known to lead to formation of binary stars, 
or, dependent on the concrete realization, even to multiple 
systems (e. g. Arreaga- Gar cia et al.||2007| . Here, the initial 
density is altered by a non-axisymmetric m = 2 perturba- 
tion in density, 



p = Po [1+Acos(20)]. 



(23) 



where the <j> is the azimuthal angle with respect to the rota- 
tion axis. Realizations of non-uniform density distributions 
for regular spaced particles, however, are not easy to achieve. 
In this case, one possible way would be the use of parti- 
cles with unequal masses. But this can lead to undesirable 
side-effects (e. g. Rosswog 2009), so this approach is not 



used by many researchers, except when the system under 



consideration is strongly centrally condensed (e. g. Arreaga- 
Garcia et al.||2010 ). For small perturbations, it is sufficient 
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M/<t> = 4 




Figure 7. Two panels show magnetic cushioning in the M/<E> = 4 run. This example shows the system at t = 1.35%, thus corresponding 
to panel 5 in Fig. |6j. The panel on the left-hand side shows column density and integrated magnetic field vectors which, for better 
visibility, were set to equal magnitude. The right-hand side shows the integrated magnetic pressure. It can be seen, that the magnetic 
field forms a cushion, which prevents the two objects from merging. This figure might be compared to the similar Fig. 11 in Price Sz 
Bate (2007), which shows the same quantities. 



to slightly perturb the initial positions of the particles, in 
order to match the desired density distribution. To do so, 
we consider the linearized continuity equation, which reads 



dp + poV • Sr = 0. 



(24) 



From this it follows, by considering spherical polar coordi- 
nates and performing a straightforward integration [using 
eq. (23)], that the perturbation in the azimuthal angle is 
given by 



Asin(20 o ) 



(25) 



For the calculations presented in this work, A = 0.1 has 
been chosen to make comparisons to other work possible. 

The cloud has a solid-body rotation with an angular 
velocity of O = 1.006 x 10~ 12 rads _1 and we fix the initial 
temperature to T = 8.4 K and the mean molecular weight 
to U- m oi = 2. With these choices, we get a sound speed of 
c s = 0.19kms _1 and initial energy ratios 



C^therm — 



-E'therm 
-Egrav | 
Erot 



= 0.26 



0.16. 



(26) 



An established quantitative measure for magnetic field sup- 
port in self-graviting fluids is the mass-to-flux ratio. For a 
spherical cloud, it is given by 



M 



7vR 2 B 



(27) 



and the critical value for this geometry (see Mouschovias & 
Spitzer 1976) in cgs units, 



: 0.125 



486 gem" VG" 



(28) 



Below this value, the magnetic field is able to support the 
cloud against gravity, while above it, gravity will dominate 
the magnetic field. At this point, we would like to stress the 
fact, that the mass-to-flux ratio is the main parameter that 
is changed in this work. Using the mass-to-flux ratio in units 
of the critical value, the magnetic field strength is given by 



B = 814 u-G 



MY 



M 



R 



4 x 10 16 cm 



(29) 



respectively. In our calculations, the initial magnetic field is 
aligned with the rotation (z) axis. 



4 RESULTS 

We chose models with several initial mass-to-flux ratios also 



considered in Price & Bate (2007), summarized in Table 
|l]), together with the corresponding field strengths B . At 
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Figure 8. Simulation results with magnetic field parallel to the rotation (z) axis, with an initial field strength of Bo = 407|U.G. This 
corresponds to a mass-to-flux ratio of M/<I> = 2, in multiples of the critical value. As previously, in horizontal direction, ^-coordinates are 
plotted, and y-coordinates are plotted vertically. Shown is the column density, here integrated along the z-axis, in physical coordinates. 
The panels display a sequence of increasing time, measured in in units of free-fall times, % = 2.4 x 10 4 yr. 



this point, we would like to emphasize, that the mass-to- 
flux ratio is given in units of the critical value throughout 
the subsequent sections. Furthermore, the table includes two 
additional columns containing the ratio of gas pressure to 
magnetic pressure, 



/3pla 

and the Alfven speeds 



Bg/(87T) 

Bo 



V&rpo 



(30) 



(31) 



respectively. 



4.1 Column density evolution 

Figure |2]), where the hydrodynamical case is shown, dis- 
plays a sequence of nine plots showing the column density 
integrated along the ^-direction, i. e. parallel to the rotation 
axis. The timing is very similar compared to Price & Bate 



Table 1. Initial parameters characterizing the magnetic field 
within the spherical cloud for each considered model. 





Bo [M-G] 


/^plasma 


va [km/s] 


oo 





OO 





20 


40.7 


39 


0.04 


10 


81.3 


9.80 


0.08 


7.5 


108.5 


5.50 


0.11 


4 


203 


1.57 


0.21 


2 


407 


0.39 


0.42 


1 


814 


0.09 


0.84 


M/$ 


is the mass 


-to-flux ratio measured 



the initial magnetic field strength, /3 p i asma is the ratio of the gas 
pressure to the magnetic pressure, and is the Alfven speed. 



2007) 



with a delay in our calculations estimated to be be- 
low 5%. Since such a disagreement in comparisons to other 
work is frequently reported in the literature, e. g. in |Com-1 



mercon et al. (2008), we regard the difference here as being 
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Figure 9. Panel of histogram plots showing the B — p relation for each initial mass-to-flux ratio where only the last snapshot prior to 
sink particle formation is considered. The grey-scaled 2D histogram indicates the value of magnetic field strength B as function of p, 
where the colour intensity is proportional to the number of particle counts within the bin considered. The solid black line is the moving 
average calculated from the contributions by the particles. The red line shows a fit assuming the relation B cx p K , where the value of k 
obtained by the fit is displayed on the lower right of each plot. Only particles with densities above po were considered in these plots. 



acceptable. Furthermore, we find that orientation and size of 
the spiral pattern agree quite well with those found in |Price| 
|&; Bate| ( [2007| ). Especially their result shown for t = 1.30t ff 
seems to coincide with our result at t = 1.35£ff. Also we no- 
tice, that rotational symmetry is well preserved. From these 
observations we conclude, that our hydrodynamical simula- 



tion results are well in agreement with other results obtained 
by other authors using a similar setup. 

The case with M/<f> — 20, which corresponds to a initial 
field strength of Bo — 40.7 |^G, is shown in Figure This is 
a comparably weak field, and so there are almost no changes 
in the spiral patterns or in timing, compared to the pure 
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hydro dynamical case. Also the onset of star formation is not 
significantly hindered. Only in the later stages, displayed in 
the lower three panels of Figure |3| , we notice a small speed- 
up in the dynamics. 

The Figures Q and ([5]), showing the case with Bo — 
81.3 nG (M/$ = 10) and B = 108.5 [iG (M/$ = 7.5), 
respectively, are more interesting. Here we see, that the col- 
lapse leads to the formation of protostars already in earlier 
stages of the simulation than in the cases considered before. 
We estimate the speed-up to be about t = 0.02 1&. Addi- 
tionally, we see that at t — 1.38 ts in the M/<f> — 10 and at 
t 1.37 ts in the M/$ 7.5 case, a third star is formed 
which in the further evolution turns out to be gravitationally 
bound to one of the other protostars, respectively. 

In order to explain the formation of a triple system in 
the latter two cases, we performed additional simulations 
in these cases. These simulations, however, do not include 
magnetic tension forces. That is, the B k B l terms in the mag- 



4.2 Relationship between B and p 



netic stress tensor, eq. (11), were neglected. In these cases 
(not shown) we find that binary systems are formed, but 
no triple system. Thus, magnetic tension could be a pos- 
sible reason for this behaviour. However, fragmentation in 
the isothermal regime is extremely sensitive to small pertur- 
bations, so small differences in the fragmentation behaviour 
should be viewed with caution. It is also known that use 
of a barotropic equation of state can severely overestimate 
the amount of secondary fragmentation compared to more 
realistic calculations where radiative heating of the gas is 



explicitly accounted fo r (e. g. lOffner et al.|[2009 Price fe| 
Sate|[2009l [Peters et al.||2010|a|b| ). 



More interesting, with respect to magnetic tension, is 
the case with mass-to-flux ratio of M/<& = 4, corresponding 
to Bo = 203 \±G, where we see a different picture (Figure |6|). 
While star formation sets in still earlier than in the previ- 
ous case, it is accompanied with a pronounced filamentary 
structure, already visible at t = 1.28 £ff. Contrary to the two 
cases considered before, only a binary system is formed here, 
where Price & Bate ( |2007| find just a single star in this case. 
This can be attributed to the 'magnetic cushioning' effect, 
shown in detail in Fig. Q at t = 1.28 tff thus corresponding 
to panel 5 in Figure This 'magnetic cushion', due to ten- 
sion forces, prevents the two protostars from merging into 
a single one, what otherwise is likely to happen given the 
close encounter of the two objects. It must be noted, how- 



ever, that Price & Bate (2007) find this 'magnetic cushion' 



only in cases where the initial magnetic field was perpen- 
dicular to the rotations axis. A possible reason for this is 
an underestimation of this effect in their calculations with 
the magnetic field parallel to the rotation axis, due to the 
restrictions of the Euler potentials in capturing certain field 
geometries. 

Figure (|8| shows the case with a very strong field, B — 
407 (M/<& = 2). Here we see a bar structure forming, 
condensing to a single star forming its central region. This 
protostar, however, is formed comparably late compared to 
the cases with weaker field strengths. Finally, we note that 
in the critical case with M/<& = 1, no star is formed at all 
(not shown). 



Observations (e. g. |Crutcher||1999| |Heiles fe Crutcher|2 005 ) 
as well as theoretical investigations and simulations (e. g. 
Mouschovias||1976| [19911 |Fiedler fe Mouschovias||1993| |DeT 



sch Mouschovias||2001| |Li et al.||2004| ) suggest a scal- 
ing behaviour of B with the density p, which is usually 
parametrized as B oc p K . 

For an isothermal (or sufficiently cooled) core, this rela- 
tion with k — 1/2 can be motivated by assuming, that mag- 
netic fields can not provide support against gravity along 
the field lines (i. e. parallel to the z axis), leading to a disc 
like morphology of the cloud in later stages of the collapse. 
By further taking the validity of flux-fre ezing in ideal MHD 



into account, one yields B oc \fpT (e. g. Heiles & Crutcher 



2005), reducing to B oc ^fp in the isothermal case. 



On the other hand, for gravity exceeding both magnetic 
and turbulent support, as well as for negligible angular ve- 



locity, the cloud could also collapse more rapidly (e. g. |Heiles 
\h Crutcher|2005|. For the ca se of a weak field and a spher- 
ical cloud, iMestel & Spitzer (1956) performed an analysis 



showing that the morphology is almost unaffected in such a 
case; furthermore they obtained a value of k — 2/3. 

In order to examine which one of these cases is realized 
in our simulation results, we show in Fig. (JqJ the depen- 
dence of the magnetic field strength B on the cloud density 
p for each of the considered cases at a time late in the evo- 
lution of the cloud, before the first sink particle is created. 
The grey-scaled two-dimensional histogram shows the mag- 
netic field strength, with the colour intensity proportional to 
the number of particles within each bin. For the calculation 
of each histogram, 200 x 200 bins were considered, respec- 
tively. The black solid curve indicates the moving average 
while the red solid curve was obtained from a fit. Our fit 
shows a power-law behaviour with a value for k which is, 
in each case, below but close to a value of 1/2 indicating 
the emergence of a disc like geometry during the collapse, in 
agreement with other studies ( Fiedler &; Mouschovias| 1993 



|Desch fc Mouschovias||2001| |Li et al.||2004| ). But note, that 
the works by Banerjee & Pudritz (2006) and Price & Bate 
(2007J report values of k 



0.6 more closely to a value of 
k = 2/3, thus indicating a more spherical collapse. However, 
our initial angular velocity Q is rather high compared to the 
latter works, who analysed the B — p relation only in an 
unperturbed setup with an gular velocities of 1.89 x 10 -13 



(iBanerjee & Pudritz||2006 ) and 1.77 x 10~ 13 rads _1 (Price 



& B ate|2007| , respectively. So in our case, flattening of the 
cloud during the collapse is expected to be enhanced com- 
pared to the latter works, especially since the effects of mag- 
netic braking, see next sub-section, are found to be rather 
weak. 



4.3 Angular momentum transport 

Also of interest is the influence of magnetic fields on the 
rotation of the cloud, usually attributed to a process called 
magnetic braking. The usual qualitative picture describing 
this process is, that torsional Alfven waves are launched into 
the ambient medium, if the latter has a different rotation 
than the cloud. Thus, the cloud is slowed down by trans- 



port of angular momentum outwards (e. g. Mouschovias & 
|Paleologou1[l979l fl980| |Mestel Par?sl|1984| |Mouschovias 
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1991). Since the ambient medium has no initial velocity at 



all, it can be expected that magnetic braking takes place in 
our models and shows measurable effects within the simu- 
lation time. To quantify the effect of magnetic braking, we 
follow the evolution of the normalized angular momentum, 
|L|/|Lo|, within the initial cloud radius R, shown in Fig. 
(10}. 

For a more quantitative analysis, we consider the 
timescale characteristic to magnetic braking, Tb, which is 
the time the outward propagating Alfven waves need to en- 
fold a fraction of the ambient gas corresponding to a mo- 
ment of inertia equal to the cloud. In the case of a spherical 
cloud, threaded by a uniform magnetic field parallel to the 



rotation axis, the braking time can be estimated by ( |McKee 
et al.||1993 ) 
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(32) 



where the primes denote values in the ambient medium. But 
note that the applicability of this classical analysis (e. g. 
Mouschovias 1991) was criticised by Hennebelle & Ciardi 
(2009), since the magnetic braking within a collapsing core 



might not be fully captured by this analysis. However, since 
we do not concentrate on the details of disc formation, but 
on the whole cloud, we expect that this analysis still gives a 
rough approximation of timescales related to magnetic brak- 
ing. By inserting numerical values, we see that the brak- 
ing times are distributed, monotonically decreasing with de- 
creasing mass-to-fiux ratio, in a range from Tb ~ 35£ff for 
M/$ = 20, to r b « 3.5 t s for M/<3> = 2. This approximation 
seems to be in good agreement with our results as displayed 
in Fig. ( |10| ) which show a rather slow braking which can 
be expected for supercritical clouds with a high initial den- 
sity ratio po/p ( McKee et al.|l993| . Furthermore, we would 
like to emphasize that our models are based on an ideal, 
one fluid MHD formulation which represents a fully ionized 
plasma, and therefore does not allow for ambipolar diffu- 
sion which limits the efficiency of magnetic braking consid- 
erably, as was pointed out already by Hosking & Whitworth 
(2004). Also of importance for the efficiency of magnetic 



braking is the initial field geometry. It was shown by |Price] 



& Bate (2007), that an initial magnetic field perpendicular 



to the rotation axis increases the efficiency of magnetic brak- 
ing substantially which is attributable to magnetic tension. 
However, || Hennebelle &; Ciardi| ( |2009| , who investigated col- 
lapse problems with magnetic fields inclined to the rotation 
axis systematically, propose that increasing the inclination 
angle of the magnetic field reduces the efficiency of magnetic 
braking. 



4.4 Numerical stability of the SPMHD algorithms 

A further aspect of our work considers the influence of di- 
vergence of the magnetic field on our results. Since non- 
vanishing divergence of the magnetic field used to impose 
serious constraints on the usability of MHD in SPH, in par- 
ticular in a star formation context, we are going to discuss 
the implications of this issue in some detail. Especially, since 
the meaning of the divergence is probably mistaken, we are 
going to clarify matters based on a similar chain of argu- 



Figure 10. Angular momentum evolution for each cloud, dis- 
tinguished by its initial mass-to-flux ratio, respectively. Shown is 
the time evolution (in %) of the normalized angular momentum, 
|L|/|Lo|, measured within the initial cloud radius R. The plot 
quantifies the decay of the initial angular momentum, dependent 
on the initial mass-to-flux ratio. 



The SPH estimator for the divergence of the magnetic 
field at the position of particle i given by, e. g., 



(V-B) 



Pi , 



■ ViW i3 



(33) 



calculates the weighted contribution of the differences of the 
magnetic field due to the N neighbouring particles to parti- 
cle i within a smoothing length h. So the magnitude of the 
divergence calculated this way, is essentially based on the 
(irregular) distribution of the particles and thus a measure 
of sub-smoothing- length fluctuations of the field. It must be 
noted, however, that this numerical divergence, which is not 
a physical divergence caused by magnetic monopoles, is also 
present in Euler potentials. The latter are free of physical 
divergence by definition, but show a mean numerical diver- 
gence, measured by the expression (/i|V • B\) / that 
can approach values of order unity during a simulation as 
reported in |Kotarba et~aT1 ( |2QQ9[ |2010| ). Additionally, the 
correction techniques employed within the present SPMHD 
scheme ensure, that the magnetic field evolution is not af- 
fected directly by the measured numerical divergence. How- 
ever, it is of course advisable to keep it as low as possible, 
and to keep track of the value of numerical divergence within 
simulation time, to ensure that irregularities in the results 
are not correlated with high values of numerical divergence. 

To quantify the analysis of the effects of numerical di- 
vergence, we plotted it in Fig. ( 11 ) as function of \B\ in a 2D 



ments as given in |Kotarba et al. (2010) 



histogram, where the intensity is proportional to the number 
of particles within a bin, with 200 x 200 bins in total. In this 
plot, the last snapshot before sink formation is considered, 
respectively, and only particles with densities larger than po 
are taken into account. It can be seen, that the numerical 
divergence is distributed on a wide range of values, as in a 
similar analysis carried out by |Kotarba et al^] |2010) for their 
systems. Additionally, no strong dependence of the (mean) 
divergence on the strength of the magnetic field, and thus 
on the density, is visible. So since the curve of the mean di- 
vergence has a trend which is qualitatively the same for all 
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Figure 11. Numerical divergence plotted as function of the magnetic field strength. The grey-scaled 2D histogram shows the values of 
/i|V • .B|/|.B| as a function of the magnetic field for particles with a density larger than po- The black solid line represents the moving 
average of the numerical divergence obtained from the individual values of the particles. As in Fig |9j, the last snapshot prior to sink 
particle formation is considered. 



initial mass-to-flux ratios, but the actual physical behaviour 
in the evolution of each setup shows huge differences, as il- 
lustrated above, we conclude that our results are meaningful 
and are not correlated to the value of numerical divergence. 



5 DISCUSSION 

We have performed a study on the influence of magnetic 
fields on collapse and fragmentation of a rotating molecular 
cloud core with an initial m — 2 density perturbation and 
the initial magnetic field aligned with the rotation axis. The 
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amplitude in each case was chosen to be 10 per cent, as 
commonly used in the literature. 

Since our approach is based on a induction equation 
formulation of SPMHD, in contrast to the Euler equations 



based approach used by Price & Bate (2007) for their star 
formation calculations, we would like to emphasize that it is 
not a priori clear that our approach should work at all for 
collapse problems, given the considerable amount of failed 
attempts in using this method (Daniel Price, private com- 
munication). Most of those attempts showed a disruptive 
behaviour at large densities, attributed to high values of nu- 
merical divergence of the magnetic field. However, since the 
details of those simulations are not known to us, we can 
only speculate what the reasons for these differences might 
be. First of all, the used parameters controlling artificial vis- 
cosity and resistivity, respectively, have values which are not 
uncommon in the literature (e. g. |Price fc Monaghan||2005| 
Dolag & Stasyszyn 2009), and thus are unlikely to have dra- 
matic effects on the global evolution of our simulations. We 
also investigated the influence of replacing high-density re- 
gions by accreting sink particles. Therefore, we performed 
simulations without any sink particles and thus followed the 
evolution of our systems as long as the global timestep al- 
lowed this at reasonable computational cost, and we have 
not recognized any signs of disruptive behaviour there. How- 
ever, since we do not know of other work that has used the 



regularization method by B0rve et al. (2001) for collapse 



simulations before, we suppose that this method, in combi- 
nation with artificial resistivity, could be more effective than 
other methods in preventing numerical divergence from cor- 
rupting the magnetic field evolution. 

Considering the column density evolution in our models, 
our results suggest an overall agreement with the findings 



from Price & Bate (2007). For very weak field strengths, as 



expected, there are almost no deviations from the pure hy- 
drodynamical case. On the other hand, for a very strong field 
with M/<f> = 2 star formation is substantially delayed and 
just a single star is formed. However, for intermediate field 
strengths with M/<3> = 10, 7.5, we actually see the formation 
of a triple system which is not present in the work by |Price] 
& Bate (2007), at least at these early stages. Using simula- 



tion without magnetic tension, the third protostar did not 
form but most probably due to small perturbations, to which 
a barotropic equations of state is very sensitive, changing 
the sub-fragmentation pattern in the considered systems. So 
we conclude, that these differences in sub-fragmentation are 
probably not very meaningful. A larger difference between 
our work and [Price &; Bate (2007) can be seen in the case 
of M/<f> — 4, where a single star is formed in their calcula- 
tions, but a binary system in our case. This difference is due 
to the 'magnetic cushioning effect' which is probably un- 
derestimated in their calculations with initial magnetic field 
parallel to the rotation axis, due to the intrinsical limitations 
of the Euler potentials. 

However, it is also instructive to compare our results 
to other findings in the literature. |Hosking &: W hitworth 
( 2004 ) investigated collapse and fragmentation using a two- 
fluid model, allowing them to model effects from non-ideal 
MHD. However, they started from sub-critical cores that be- 
came critical during the evolution via ambipolar diffusion. 
So it must be noted, that their investigations are quite dif- 
ferent from our approach. They found no fragmentation in 



their magnetized models, but due to the different approaches 
used, it is difficult to relate their findings to our results. 

Furthermore, we would like to mention the work by 
Machida et al. (2005b), who investigated a large range of 



parameters in their fragmentation problems. However, since 
they started out from a filament with more complex initial 
perturbations in the density as well as in the magnetic field 
itself, their initial conditions are very different from those 
used in this work. So any comparison can be only of quali- 
tative nature. Their characterizing parameter uj corresponds 
to \//?rot = 0.4, while their a is equal to v\/cl. Thus, our 
models in a range from M/<£ = 20 to M/<& = 4 are located 
in the 'vertical collapse region' in their Figure 10, were frag- 
mentation is possible according to their analysis. The very 
strong field models with M/4> = 2 to M/4> = 1 are outside 
the horizontal range of their Figure 10, but the rather ex- 
treme values of a = 4.88 and a = 19.54, respectively, lead us 
to the speculation that they would be located in the region 
were no fragmentation occurs. Therefore, we find that our 
results globally agree with those found by Machida et al. 
(12005b). 



A comparison to Hennebelle & Teyssier ( 2008 ) is quite 



difficult, since they use a different barotropic equation of 
state with a critical density of p cr it = 10 _13 g cm -3 , the latter 
being one order of magnitude higher than ours. Additionally, 
their f3 vo t has a value of 0.045 lower than in our models. In 
their weak perturbation models using A = 0.1, they find no 
fragmentation for M/4> ^ 20, thus all of their models form 
a single star. This is, with the exception of the M/<& = 2 
case, in disagreement with our findings. However, this not 
surprising because of the differences to their initial setup. 

Ziegler (2005) and Fromang et al. (2006), investigated 
the case M/$ = 2 using the same initial conditions, which 



are very similar to those used by Hennebel le &; Teyssier 
2008|). They also u sed the same equation of state as [Hen- 
nebelle & Teyssierl (120081) and f3 ro t = 0.045, thus only a 



qualitative comparison is possible to our work. Ziegler ( 2005 ) 
finds formation of a binary in this case, while Fromang et al.| 
( 2006 ) get different results depending on the flux solver used, 
namely no binary with the Lax-Friedrich solver and a binary 
with the Roe solver. However, the latter binary merged to a 
single fragment shortly thereafter. Our results for M/<& = 2 
show no sign of binary formation, so considering this partic- 
ular case our results show more similarities with Fromang 
et al.| ( |2Q06| ) than with |Ziegler| ( |2005| ). 

Furthermore, we also investigated for each mass-to- flux 
ration the dependence of the magnetic field strength B on 
the density p within the final stages of collapse within a 
core. Our results, showing a value of k close to 1/2 in the 
power-law relation B oc p K , are well in agreement with 
a picture with vanishing magnetic support parallel to the 
symmetry axis. Thus, the cloud finally ends in a disc-like 
morphology, independent of the initial mass-to- flux ratio. 
Such a behaviour is also frequently reported in the litera- 



ture (|Mouschovias|1976l|1991||Fiedler fe Mouschovias|199"3 



Desch fe Mouschov ias 2001; " Li et al.|2004 ) 

An additional investigation concerned the angular mo- 
mentum transport, yielding that magnetic braking is weak 
in the models we considered, as could be expected from an- 



alytical reasoning dMouschovias fe Pa leologou 197 9] |1980 



Mestel fe Paris|1984||Mouschovias|1991||McKee et al.|1993f, 
and from non-ideal MHD simulations carried out by Hosking 
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& Whitworth ( 2004 ) . The influence of the initial field geom- 



etry on the efficiency of magnetic braking is currently under 
discussion. iPrice &; Bate (2007) advocate an increased effi- 



ciency with an an initial field perpendicular to the rotation 



axis, but Hennebelle & Ciardi (2009) propose the opposite 



However, we would like to emphasize, that we consider the 
whole cloud in this analysis and pay no attention on the im- 
pact of magnetic braking on disc formation. Thus we regard 
the recent criticism of the classical analysis by Hennebelle 
& Ciardi (2009) as not influential to our analysis. 



Finally, we would like to stress the fact, that according 
to our analysis of numerical divergence, we do expect that 
our results are not corrupted by artefacts and that therefore 
our results show the correct physical behaviour within our 
systems. 
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6 SUMMARY 

In this work, we carried out magnetohydrodynamical com- 
puter simulations of the collapse of molecular cloud cores, 
initially disturbed with m — 2 density perturbations. 

The method is based on a formulation of smoothed 
particle magnetohydrodynamics (SPMHD) that evolves the 
magnetic field directly via the induction equation and thus 
does not utilize any form of scalar or vector potentials. Sta- 
bility and noise reduction are ensured by techniques imple- 
mented by Dolag & Stasyszyn (2009), which have, to the 



best of our knowledge, not yet been applied in this combi- 
nation to star formation problems. 

From the results of this work, we draw several main con- 
clusions. First, we find that our formulation of SPMHD did 
well in reproducing essential features obtained in other work 
with similar initial conditions, but using different methods. 
Thus we conclude that our approach is a viable scheme to at- 
tack star formation problems. Second, our results show good 



global agreement with the work by Price & Bate (2007), 



with the exception of cases with higher field strength where 
magnetic tensions aids binary fragmentation via the 'mag- 
netic cushioning effect' in our simulations. This effect is not 



present in the corresponding results in Price & Bate (2007), 



most probably due to limitations of the Euler potentials ap- 
proach in representing certain geometries of the magnetic 
field. 
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